# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <math.h>
# include <stdlib.h>

# ifdef _OMP_VEL_ADVANCE_
# include <omp.h>
# endif

extern long int solve_poisson(struct grid *g, double ****rhs, unsigned short iRK);

extern void vel_bc(struct grid *g, short tflag, unsigned short iRK);
extern void press_bc(struct grid *g, short tflag, unsigned short iRK);

extern double dens(double T);
extern double kvisc(double T);

static double getA(int i, int j, int k, struct grid *g, unsigned short isec);
static double getB(int i, int j, int k, struct grid *g, unsigned short isec);
static double getC(int i, int j, int k, struct grid *g, unsigned short isec);

struct diverg_max
     {
      short im, jm, km;
      double divmax;
     };

static void divergence(struct grid *g, double ***diverg, struct diverg_max *mdiv, unsigned short isec);

extern unsigned int vel_advance(struct grid *g, unsigned int *pcnt_max, struct arr_struct *depVar, struct arr_struct (*memFace)[3])
{
 unsigned int pcnt, pcnt_subttl;

 short i, j, k, _i, _j, _k, nr, nth, nz, isec;
 double density0, dt, *r, *th, *z, ***d[Nsec], ***D0[Nsec],\
        ***_u[Nsec][nRK+1], ***_v[Nsec][nRK+1], ***_w[Nsec][nRK+1], ***_p[Nsec][nRK],\
        ***p, ***u0, ***v0, ***w0, ***p0, ***u1, ***v1, ***w1, ***p1, ***u2, ***v2, ***w2, ***p2, ***u3, ***v3, ***w3,\
        ***_A[Nsec][nRK], ***_B[Nsec][nRK], ***_C[Nsec][nRK], ***A0, ***B0, ***C0, ***A1, ***B1, ***C1, ***A2, ***B2, ***C2,\
        ***divH[Nsec][nRK-1], ***divH0, ***divH1, divH2_ijk, ***div2p0[Nsec], div2p1_ijk;

// struct timeb t1, t2;

 th = g->th; dt = g->dt; nth = g->nth;

 *pcnt_max = pcnt_subttl = 0;
 density0 = g->density0;

/* Memory for all section of the computational domain. */
  for(isec = 0; isec < Nsec; isec++)
    {
     _u[isec][0] = g->u[isec];         _v[isec][0] = g->v[isec];         _w[isec][0] = g->w[isec];          _p[isec][0] = g->p[isec];
     _u[isec][1] = depVar[isec].q[4];  _v[isec][1] = depVar[isec].q[5];  _w[isec][1] = depVar[isec].q[6];   _p[isec][1] = depVar[isec].q[7];
     _u[isec][2] = depVar[isec].q[8];  _v[isec][2] = depVar[isec].q[9];  _w[isec][2] = depVar[isec].q[10];  _p[isec][2] = depVar[isec].q[11];
     _u[isec][3] = depVar[isec].q[12]; _v[isec][3] = depVar[isec].q[13]; _w[isec][3] = depVar[isec].q[14];

     D0[isec] = depVar[isec].q[11]; divH[isec][0] = depVar[isec].q[12]; divH[isec][1] = depVar[isec].q[13]; div2p0[isec] = depVar[isec].q[14];

     for(i = 0; i < nRK; i++)
         {
          _A[isec][i] = memFace[isec][0].q[i];
          _B[isec][i] = memFace[isec][1].q[i];
          _C[isec][i] = memFace[isec][2].q[i];
         }
    }

/* Calculating divergence field at the current n-th time-step. It gets filled in D0[isec][][][]. */
  for(isec = 0; isec < Nsec; isec++) divergence(g, D0[isec], (struct diverg_max *)0, isec);

/* --------------------------- Calculation of u(1), v(1), w(1), and p(1) ---------------------------------- */
  for(isec = 0; isec < Nsec; isec++)
      {
       nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
       A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0];
       u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0]; p0 = _p[isec][0];
       u1 = _u[isec][1]; v1 = _v[isec][1]; w1 = _w[isec][1];

#      ifdef _OMP_VEL_ADVANCE_
       #pragma omp parallel private(i,j,k)
#      endif
          {
#          ifdef _OMP_VEL_ADVANCE_
            #pragma omp for schedule(dynamic,1)
#          endif
           for(k = 0; k <= nz; k++) for(i = 0; i <= nr; i++) for(j = 0; j <= nth; j++)
              {
               if(j < nth && k < nz) A0[i][j][k] = getA(NrG+i-1, NthG+j, NzG+k, g, isec);
               if(i < nr && k < nz) B0[i][j][k] = getB(NrG+i, NthG+j-1, NzG+k, g, isec);
               if(i < nr && j < nth) C0[i][j][k] = getC(NrG+i, NthG+j, NzG+k-1, g, isec);
              }

#          ifdef _OMP_VEL_ADVANCE_
            #pragma omp for schedule(dynamic,1)
#          endif
           for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
              {
               u1[i][j][k] = u0[i][j][k]\
                 + c2*dt*( A0[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p0[i+1][j][k] - p0[i][j][k])/(r[i+2]-r[i]) );

               v1[i][j][k] = v0[i][j][k]\
                 + c2*dt*( B0[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p0[i][j+1][k] - p0[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) );

               w1[i][j][k] = w0[i][j][k]\
                 + c2*dt*( C0[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p0[i][j][k+1] - p0[i][j][k])/(z[k+2]-z[k]) );

               if(i == NrG+nr-1) u1[i][j][k] = 0;
               if(isec == 1 && k == NzG+nz-1) w1[i][j][k] = 0;
              }
          }

       g->u[isec] = u1; g->v[isec] = v1; g->w[isec] = w1;
      }

  vel_bc(g, 0, 1);

/* _u[1], _v[1], _w[1] is used in "press_bc()". Therefore, the memory used by these can't be distroyed. */
  for(isec = 0; isec < Nsec; isec++)
      {
       nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
       A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0]; p0 = _p[isec][0]; divH0 = divH[isec][0];
       A1 = _A[isec][1]; B1 = _B[isec][1]; C1 = _C[isec][1]; p1 = _p[isec][1]; divH1 = divH[isec][1];

       d[isec] = depVar[isec].q[8]; p = _p[isec][0];

#      ifdef _OMP_VEL_ADVANCE_
       #pragma omp parallel private(i,j,k,_i,_j,_k)
#      endif
         {
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic,1)
#         endif
          for(k = 0; k <= nz; k++) for(i = 0; i <= nr; i++) for(j = 0; j <= nth; j++)
              {
               if(j < nth && k < nz) A1[i][j][k] = getA(NrG+i-1, NthG+j, NzG+k, g, isec);
               if(i < nr && k < nz)  B1[i][j][k] = getB(NrG+i, NthG+j-1, NzG+k, g, isec);
               if(i < nr && j < nth) C1[i][j][k] = getC(NrG+i, NthG+j, NzG+k-1, g, isec);
              }

       /* Calculating RHS of the PPE-1. */
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic,1) nowait
#         endif
          for(_k = 0; _k < nz; _k++) for(k = _k+NzG, _i = 0; _i < nr; _i++) for(i = _i+NrG, _j = 0, j = NthG; _j < nth; j++, _j++)
               {
                div2p0[isec][i][j][k] = DIV2P_ijk;

                divH0[i][j][k] = 2*(( r[i+1]*A0[_i+1][_j][_k] - ( ((int)(r[i]*10000) == 0) ? A0[_i][_j][_k]:r[i]*A0[_i][_j][_k]))/(r[i+1]-r[i])\
                           + (B0[_i][_j+1][_k] - B0[_i][_j][_k])/(th[j+1]-th[j]) )/(r[i]+r[i+1])\
                           + (C0[_i][_j][_k+1] - C0[_i][_j][_k])/(z[k+1]-z[k]);

                divH1[i][j][k] = 2*( (r[i+1]*A1[_i+1][_j][_k] - (((int)(r[i]*10000) == 0) ? A1[_i][_j][_k]:r[i]*A1[_i][_j][_k]))/(r[i+1]-r[i])\
                          + (B1[_i][_j+1][_k] - B1[_i][_j][_k])/(th[j+1]-th[j]) )/(r[i]+r[i+1])\
                          + (C1[_i][_j][_k+1] - C1[_i][_j][_k])/(z[k+1]-z[k]);

                d[isec][i][j][k] = density0*( divH1[i][j][k]\
                                            + (D0[isec][i][j][k]/dt + a31*(divH0[i][j][k] - div2p0[isec][i][j][k]/density0))/a32 );
               }

       /* Copying _p[0] into _p[1] to provide it as initial condition for pressure while solving the PPE. */
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic,2)
#         endif
          for(k = 0; k < 2*NzG+nz; k++) for(i = 0; i < 2*NrG+nr; i++) for(j = 0; j < 2*NthG+nth; j++)
             p1[i][j][k] = p0[i][j][k];
         }
       g->p[isec] = p1; /* Used in solve_poisson(). */
      }


/* Solving the PPE(1) to calculate p(1). */
   if(g->tptr != (FILE *)0) fprintf(g->tptr, "(PPE - 1)\n");
   pcnt = solve_poisson(g, d, 1);

   pcnt_subttl += pcnt;
   if(pcnt > *pcnt_max) *pcnt_max = pcnt;
/* ------------------------------------------------------------------------------------------------------------------------- */


/* ------------------------------------- Calculation of u(2), v(2), w(2), and p(2) ----------------------------------------- */
/* Calculating u(2), v(2), and w(2). */

  for(isec = 0; isec < Nsec; isec++)
     {
      nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec]; p0 = _p[isec][0];  p1 = _p[isec][1];
      u2 = _u[isec][2]; v2 = _v[isec][2]; w2 = _w[isec][2]; u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0];
      A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0]; A1 = _A[isec][1]; B1 = _B[isec][1]; C1 = _C[isec][1];

#     ifdef _OMP_VEL_ADVANCE_
      #pragma omp parallel for private(i,j,k) schedule(dynamic,2)
#     endif
      for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
          {
           u2[i][j][k] = u0[i][j][k]\
                + a31*dt*( A0[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p0[i+1][j][k] - p0[i][j][k])/(r[i+2] - r[i]) )\
                + a32*dt*( A1[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p1[i+1][j][k] - p1[i][j][k])/(r[i+2] - r[i]) );

           v2[i][j][k] = v0[i][j][k]\
                + a31*dt*( B0[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p0[i][j+1][k] - p0[i][j][k])/( (r[i]+r[i+1])*(th[j+2]-th[j]) ) )\
                + a32*dt*( B1[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p1[i][j+1][k] - p1[i][j][k])/( (r[i]+r[i+1])*(th[j+2]-th[j]) ) );

           w2[i][j][k] = w0[i][j][k]\
                + a31*dt*( C0[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p0[i][j][k+1] - p0[i][j][k])/(z[k+2]-z[k]) )\
                + a32*dt*( C1[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p1[i][j][k+1] - p1[i][j][k])/(z[k+2]-z[k]) );

           if(i == NrG+nr-1) u2[i][j][k] = 0;
           if(isec == 1 && k == NzG+nz-1) w2[i][j][k] = 0;
          }

      g->u[isec] = u2; g->v[isec] = v2; g->w[isec] = w2;
     }


  vel_bc(g, 0, 2);


/* _u[2], _v[2], _w[2] is used in "press_bc()". Therefore, the memory used by these can't be distroyed. */

  for(isec = 0; isec < Nsec; isec++)
     {
      nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
      A2 = _A[isec][2]; B2 = _B[isec][2]; C2 = _C[isec][2]; p1 = _p[isec][1]; p2 = _p[isec][2];
      d[isec] = depVar[isec].q[4]; p = _p[isec][1]; divH0 = divH[isec][0]; divH1 = divH[isec][1];

#     ifdef _OMP_VEL_ADVANCE_
      #pragma omp parallel private(i,j,k,_i,_j,_k,div2p1_ijk,divH2_ijk)
#     endif
          {
#          ifdef _OMP_VEL_ADVANCE_
           #pragma omp for schedule(dynamic,1)
#          endif
           for(k = 0; k <= nz; k++) for(i = 0; i <= nr; i++) for(j = 0; j <= nth; j++)
               {
                if(j < nth && k < nz) A2[i][j][k] = getA(NrG+i-1, NthG+j, NzG+k, g, isec);
                if(i < nr && k < nz)  B2[i][j][k] = getB(NrG+i, NthG+j-1, NzG+k, g, isec);
                if(i < nr && j < nth) C2[i][j][k] = getC(NrG+i, NthG+j, NzG+k-1, g, isec);
               }

#          ifdef _OMP_VEL_ADVANCE_
           #pragma omp for schedule(dynamic,1)
#          endif
           for(_k = 0; _k < nz; _k++) for(k = _k+NzG, _i = 0; _i < nr; _i++) for(i = _i+NrG, _j = 0, j = NthG; _j < nth; j++, _j++)
               {
                div2p1_ijk = DIV2P_ijk;
                divH2_ijk = 2*( (r[i+1]*A2[_i+1][_j][_k] - (((int)(r[i]*10000) == 0) ? A2[_i][_j][_k]:r[i]*A2[_i][_j][_k]))/(r[i+1]-r[i])\
                     + (B2[_i][_j+1][_k] - B2[_i][_j][_k])/(th[j+1]-th[j]) )/(r[i]+r[i+1])\
                     + (C2[_i][_j][_k+1] - C2[_i][_j][_k])/(z[k+1]-z[k]);

                d[isec][i][j][k] = density0*( divH2_ijk + (D0[isec][i][j][k]/dt + a41*(divH0[i][j][k] - div2p0[isec][i][j][k]/density0)\
                                                         + a42*(divH1[i][j][k] - div2p1_ijk/density0))/a43 );
               }

        /* Copying _p[1] into _p[2] to provide it as initial condition for pressure while solving the PPE. */
#          ifdef _OMP_VEL_ADVANCE_
           #pragma omp for schedule(static)
#          endif
        /* IMPORTANT: The "D0" and "p2" are sharing same memory space i.e. the above calculation of "d" must be over before "p2"
           gets written below. i.e. we can't use the "nowait" flag in the above parallelization. */
           for(k = 0; k < 2*NzG+nz; k++) for(i = 0; i < 2*NrG+nr; i++) for(j = 0; j < 2*NthG+nth; j++)
            p2[i][j][k] = p1[i][j][k];
          }
      g->p[isec] = p2;
     }


/* Solving the PPE(2) to calculate p(2). */
   if(g->tptr != (FILE *)0) fprintf(g->tptr, "(PPE - 2)\n");
   pcnt = solve_poisson(g, d, 2);

   pcnt_subttl += pcnt;
   if(pcnt > *pcnt_max) *pcnt_max = pcnt;
/* ---------------------------------------------------------------------------------------------------------- */

/* ================================= Calculation of u, v, w, and p at next time-step ========================== */
/* Calculation of u, v, and w at next time-step. */
  for(isec = 0; isec < Nsec; isec++)
      {
       nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
       u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0];
       u3 = _u[isec][3]; v3 = _v[isec][3]; w3 = _w[isec][3];
       A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0]; p0 = _p[isec][0];
       A1 = _A[isec][1]; B1 = _B[isec][1]; C1 = _C[isec][1]; p1 = _p[isec][1];
       A2 = _A[isec][2]; B2 = _B[isec][2]; C2 = _C[isec][2]; p2 = _p[isec][2];

#      ifdef _OMP_VEL_ADVANCE_
       #pragma omp parallel for private(i,j,k) schedule(dynamic,2)
#      endif
       for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
           {
            u3[i][j][k] = u0[i][j][k]\
                  + b1*dt*( A0[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p0[i+1][j][k] - p0[i][j][k])/(r[i+2]-r[i]) )\
                  + b2*dt*( A1[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p1[i+1][j][k] - p1[i][j][k])/(r[i+2]-r[i]) )\
                  + b3*dt*( A2[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p2[i+1][j][k] - p2[i][j][k])/(r[i+2]-r[i]) );

            v3[i][j][k] = v0[i][j][k]\
                  + b1*dt*( B0[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p0[i][j+1][k] - p0[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) )\
                  + b2*dt*( B1[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p1[i][j+1][k] - p1[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) )\
                  + b3*dt*( B2[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p2[i][j+1][k] - p2[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) );

            w3[i][j][k] = w0[i][j][k]\
                  + b1*dt*( C0[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p0[i][j][k+1] - p0[i][j][k])/(z[k+2]-z[k]) )\
                  + b2*dt*( C1[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p1[i][j][k+1] - p1[i][j][k])/(z[k+2]-z[k]) )\
                  + b3*dt*( C2[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p2[i][j][k+1] - p2[i][j][k])/(z[k+2]-z[k]) );

            if(i == NrG+nr-1) u3[i][j][k] = 0;
            if(isec == 1 && k == NzG+nz-1) w3[i][j][k] = 0;
           }

       g->u[isec] = u3; g->v[isec] = v3; g->w[isec] = w3;
      }

  vel_bc(g, 1, 0);

/* Storing mean velocities and copying the newly evaluated velocities to _u[0], _v[0], and _w[0]. */
  for(isec = 0; isec < Nsec; isec++)
     {
      nr = g->nr[isec]; nz = g->nz[isec];
      p0 = _p[isec][0]; p1 = _p[isec][1]; p2 = _p[isec][2];
      u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0];
      u3 = _u[isec][3]; v3 = _v[isec][3]; w3 = _w[isec][3];

   /* Memory to store mean velocities at time = t + dt/2. */
      g->um[isec] = depVar[isec].q[4]; g->vm[isec] = depVar[isec].q[5]; g->wm[isec] = depVar[isec].q[6];

#     ifdef _OMP_VEL_ADVANCE_
      #pragma omp parallel private(i,j,k)
#     endif
         {
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(static) nowait
#         endif
          for(k = 0; k < nz+2*NzG; k++) for(i = 0; i < nr+2*NrG; i++) for(j = 0; j < nth+2*NthG; j++)
              {
            /* Storing mean-velocities. */
               g->um[isec][i][j][k] = 0.5*(u0[i][j][k] + u3[i][j][k]);
               g->vm[isec][i][j][k] = 0.5*(v0[i][j][k] + v3[i][j][k]);
               g->wm[isec][i][j][k] = 0.5*(w0[i][j][k] + w3[i][j][k]);

           /* Copying newly evaluated velocities. */
              u0[i][j][k] = u3[i][j][k];
              v0[i][j][k] = v3[i][j][k];
              w0[i][j][k] = w3[i][j][k];
             }

#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(static)
#         endif
          for(k = NzG; k < nz + NzG; k++) for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
          /* Pressure at the next time-step will be given by _p[0]. */
             p0[i][j][k] = d1*p0[i][j][k] + d2*p1[i][j][k] + d3*p2[i][j][k];
         }
      g->u[isec] = u0; g->v[isec] = v0; g->w[isec] = w0; g->p[isec] = p0;
     }

  press_bc(g, 1, 0);
/* ---------------------------------------------------------------------------------------------------------------- */

 return pcnt_subttl;

/* q[0]~q[6] are engagged. These can't be used in other functions. */
}


/* i,j,k -th CELL for isec-th section. */
static double getA(int i, int j, int k, struct grid *g, unsigned short isec)
 {
  double ***u, ***v, ***w, ***T, *r, *th, *z, kvisc0;

  u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; T = g->T[isec]; r = g->r[isec];  th = g->th;  z = g->z[isec];
  kvisc0 = g->kvisc0;

  return ((int)(r[i+1]*10000) == 0) ? rA_CENTER : AA;
 }


static double getB(int i, int j, int k, struct grid *g, unsigned short isec)
 {
  double ***u, ***v, ***w, ***T, *r, *th, *z, kvisc0;
  u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; T = g->T[isec];  r = g->r[isec];  th = g->th;  z = g->z[isec];
  kvisc0 = g->kvisc0;

  return BB;
 }


static double getC(int i, int j, int k, struct grid *g, unsigned short isec)
 {
  double ***u, ***v, ***w, ***T, *r, *th, *z, kvisc0, density0;

  u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; T = g->T[isec];  r = g->r[isec];  th = g->th;  z = g->z[isec];
  kvisc0 = g->kvisc0; density0 = g->density0;

  return CC;
 }


static void divergence(struct grid *g, double ***diverg, struct diverg_max *mdiv, unsigned short isec)
   {
    double *r, *th, *z, ***u, ***v, ***w, *div;
    short i, j, k, tid, nr, nz, nth, nTRD;
    struct diverg_max *tddiv;

    u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; r = g->r[isec]; th = g->th; z = g->z[isec];
    nr = g->nr[isec]; nth = g->nth; nz = g->nz[isec];

    nTRD = 1;
#   ifdef _OMP_VEL_ADVANCE_
    nTRD = g->nTRD;
#   endif

    tddiv = (struct diverg_max *)malloc(sizeof(struct diverg_max)*nTRD);
    div = (double *)malloc(sizeof(double)*nTRD);

    tid = 0;
#   ifdef _OMP_VEL_ADVANCE_
    #pragma omp parallel private(i,j,k,tid)
#   endif
         {
#         ifdef _OMP_VEL_ADVANCE_
          tid = omp_get_thread_num();
#         endif

          tddiv[tid].im = -1;
          tddiv[tid].jm = -1;
          tddiv[tid].km = -1;
          tddiv[tid].divmax = 0.0;

#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic, 3)
#         endif
          for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < nth+NthG; j++) for(k = NzG; k < nz+NzG; k++)
             {
              div[tid] = DIV_ijk;

           /* Filtering out the maximum divergence. */
              if((div[tid] < 0 ? -div[tid]:div[tid]) > (tddiv[tid].divmax < 0 ? -tddiv[tid].divmax:tddiv[tid].divmax))
                 {
                  tddiv[tid].divmax = div[tid];
                  tddiv[tid].im = i;
                  tddiv[tid].jm = j;
                  tddiv[tid].km = k;
                 }

              diverg[i][j][k] = div[tid];
             }
         }

    free(div);

 /* Calculating the absolute maximum value of the divergence. */
    for(i = 0, tid = 1; tid < nTRD; tid++)
    if((tddiv[tid].divmax < 0 ? -tddiv[tid].divmax:tddiv[tid].divmax) > (tddiv[i].divmax < 0 ? -tddiv[i].divmax:tddiv[i].divmax)) i = tid;

    if(g->tptr != (FILE *)0) fprintf(g->tptr, "\tdivmax(%d, %d, %d | %d) = %E\n", tddiv[i].im, tddiv[i].jm, tddiv[i].km, isec, tddiv[i].divmax);

    if( mdiv != (struct diverg_max *)0 )
      { mdiv->divmax = tddiv[i].divmax; mdiv->im = tddiv[i].im; mdiv->jm = tddiv[i].jm; mdiv->km = tddiv[i].km; }

    free(tddiv);

    return;
   }
